##########################################################################
## Replication syntax for Hegewald & Schraff (2025, CPS):
## "Is the Urban-Rural Divide Affectively Polarised? Comparative Evidence From Nine European Countries"
##########################################################################

# R version: 4.5.0 
# Platform: x86_64-apple-darwin20

# 1. Load packages #######################################################
library(here)
library(tidyverse)
library(magrittr)
library(psych)
library(ggpubr)
library(ggeffects)
library(patchwork)
library(ggrepel)
library(ggh4x)
library(stargazer)
library(modelsummary)
library(interflex)
library(reshape2)
library(kableExtra)
library(marginaleffects)
library(lme4)
library(performance)
library(lmerTest)

# Package version info: 
# here_1.0.1 tidyverse_2.0.0 magrittr_2.0.3 psych_2.5.3 ggpubr_0.6.0 ggeffects_2.2.1 patchwork_1.3.0 
# ggrepel_0.9.6 ggh4x_0.3.1 stargazer_5.2.3 modelsummary_2.3.0 interflex_1.2.6 reshape2_1.4.4 
# kableExtra_1.4.0 marginaleffects_0.27.0 lme4_1.1.37 performance_0.13.0 lmerTest_3.1.3

# 2. Basic setup #########################################################
rm(list=ls())
i_am("Hegewald_Schraff_CPS_replication.R")

# 3. Import data ##########################################################
dat_full <- read.csv((here("dat_fin.csv")), sep=",") 
CHES_ts <- read.csv((here("CHES_ts.csv")), sep=",", stringsAsFactors = TRUE)

# 4. Prepare data ##########################################################
dat_full %<>% filter(gender!="non-binary") ### 11 observations dropped

### Create dependent variables based on CHES 2023
dat_full$galtan_vote <- dat_full$galtan23 

cor(dat_full$galtan19, dat_full$galtan23,  method = "pearson", use = "complete.obs")

dat_full$radicalright <- ifelse(dat_full$family23 == "TAN", 1, 0) 
dat_full$green <- ifelse(dat_full$family23 == "green", 1, 0) 

### Create thermometer differential
dat_full$thermo.diff <- dat_full$thermo.in-dat_full$thermo.out

### Create trait-rating differential 
### trait-ratings in
cor(dat_full[,c(22:26)])
dat_full$self.in <- dat_full$self.in*-1+6
dat_full$hypo.in <- dat_full$hypo.in*-1+6
cor(dat_full[,c(22:26)])

alpha(dat_full[,c(22:26)]) ### 0.74 
dat_full$trait.in <- (dat_full$int.in + dat_full$open.in + dat_full$hon.in + dat_full$self.in + dat_full$hypo.in)/5

### trait-ratings out
cor(dat_full[,c(27:31)])
dat_full$self.out <- dat_full$self.out*-1+6
dat_full$hypo.out <- dat_full$hypo.out*-1+6
cor(dat_full[,c(27:31)])

alpha(dat_full[,c(27:31)]) ### 0.66

dat_full$trait.out <- (dat_full$int.out + dat_full$open.out + dat_full$hon.out + dat_full$self.out + dat_full$hypo.out)/5

dat_full$trait.diff <- dat_full$trait.in-dat_full$trait.out

### place-based resentment
cor(dat_full[,c(32:36)])
alpha(dat_full[,c(32:36)]) ### 0.83

dat_full$resent_munis <- (dat_full$eco + dat_full$rep1 + dat_full$rep2 + dat_full$cult1 + dat_full$cult2)/5

dat_full$resent_munis.std <- as.numeric(scale(dat_full$resent_munis)) 
summary(dat_full$resent_munis.std)
sd(dat_full$resent_munis.std)

### place-based identity
dat_full$urid_urb <- ifelse(dat_full$rur=="urban", dat_full$urid, NA)
dat_full$rurid_rur <- ifelse(dat_full$rur=="rural", dat_full$rurid, NA)

dat_full$id <- ifelse(dat_full$rur=="rural", dat_full$rurid_rur, dat_full$urid_urb)
dat_full$id.std <- as.numeric(scale(dat_full$id))

### standardize variables
dat_full$thermo.diff.std <- as.numeric(scale(dat_full$thermo.diff))
dat_full$thermo.in.std <- as.numeric(scale(dat_full$thermo.in))
dat_full$thermo.out.std <- as.numeric(scale(dat_full$thermo.out))
dat_full$trait.diff.std <- as.numeric(scale(dat_full$trait.diff))

dat_full$age.std <- as.numeric(scale(dat_full$age))
dat_full$lr.std <- as.numeric(scale(dat_full$lr))

dat_full$plid.std <- as.numeric(scale(dat_full$plid))

dat_full$immi.std <- as.numeric(scale(dat_full$immi))
dat_full$reg_unemp.std <- as.numeric(scale(dat_full$reg_unemp))
dat_full$pop_delta.std <- as.numeric(scale(dat_full$pop_delta))

### set base manually
dat_full$rur <- ifelse(dat_full$rur=="urban", "Urban residence", "Rural residence")
dat_full$rur <- as.factor(dat_full$rur)
dat_full <- within(dat_full, rur <- relevel(rur, ref = "Urban residence"))

dat_full$rur_cont <- as.factor(dat_full$rur_cont)
dat_full$rur_cont <- ordered(dat_full$rur_cont, levels = c("Very urban", "Rather urban", "Rather rural", "Very rural"))

dat_full$rur_cont <-  factor(dat_full$rur_cont, ordered = FALSE)
dat_full <- within(dat_full, rur_cont <- relevel(rur_cont, ref = "Very urban"))

dat_full$gender <- as.factor(dat_full$gender)
dat_full <- within(dat_full, gender <- relevel(gender, ref = "male"))

dat_full$educ <- as.factor(dat_full$educ)
dat_full <- within(dat_full, educ <- relevel(educ, ref = "Low"))

dat_full$migbck <- as.factor(dat_full$migbck)
dat_full <- within(dat_full, migbck <- relevel(migbck, ref = "No"))

### Change country labels to names 
dat_full$cntry[dat_full$cntry=="DK"] <- "Denmark"
dat_full$cntry[dat_full$cntry=="DE"] <- "Germany"
dat_full$cntry[dat_full$cntry=="EL"] <- "Greece"
dat_full$cntry[dat_full$cntry=="ES"] <- "Spain"
dat_full$cntry[dat_full$cntry=="FR"] <- "France"
dat_full$cntry[dat_full$cntry=="IT"] <- "Italy"
dat_full$cntry[dat_full$cntry=="CZ"] <- "Czech Republic"
dat_full$cntry[dat_full$cntry=="HU"] <- "Hungary"
dat_full$cntry[dat_full$cntry=="PL"] <- "Poland"
dat_full$cntry <- as.factor(dat_full$cntry)

### residential movement variable
dat_full$resbio <- factor(dat_full$resbio, levels = c("A big city", "The suburbs or outskirts of a big city", 
                                                      "A town or small city", "A country village", "A farm or home in the countryside"))

dat_full$res <- factor(dat_full$res, levels = c("A big city", "The suburbs or outskirts of a big city", 
                                                "A town or small city", "A country village", "A farm or home in the countryside"))

dat_full$resbio_num <- as.numeric(dat_full$resbio)
dat_full$res_num <- as.numeric(dat_full$res)
dat_full$movement <- dat_full$resbio_num-dat_full$res_num

dat_full$urban_move <- "No"
dat_full$urban_move[dat_full$res_num<dat_full$resbio_num] <- "Yes"
table(dat_full$urban_move)

dat_full$urban_move <- as.factor(dat_full$urban_move)
dat_full <- within(dat_full, urban_move <- relevel(urban_move, ref = "No"))

dat_full$rural_move <- "No"
dat_full$rural_move[dat_full$res_num>dat_full$resbio_num] <- "Yes"
table(dat_full$rural_move)

dat_full$rural_move <- as.factor(dat_full$rural_move)
dat_full <- within(dat_full, rural_move <- relevel(rural_move, ref = "No"))

table(dat_full$rural_move, dat_full$urban_move)

### define critical values
interval1 <- -qnorm((1-0.95)/2)  # 95% multiplier
interval2 <- -qnorm((1-0.99)/2)  # 99% multiplier

# 5. Main analysis #########################################################
# 5.1 Figure 1 #############################################################
grp.mean <- dat_full %>%
  group_by(cntry, rur) %>%
  dplyr::summarize(Mean = mean(thermo.diff, na.rm=TRUE))

png(file=here("Figure1.png"),
    width = 8, height = 6, units = 'in', res = 800)
ggplot(dat_full, aes(x = thermo.diff, fill = rur)) + 
  geom_density(alpha = .5) +
  scale_fill_manual(name = "Subsample:", values = c("black", "grey")) + 
  labs(y = "Density",
       x = "Thermometer differential") + 
  facet_wrap(vars(cntry)) +
  geom_vline(data = grp.mean, aes(xintercept = Mean, color = rur),
             linetype = "dashed") + 
  scale_color_manual(name = "Subsample:", values = c("grey", "black")) +
  theme_pubr() +
  theme(legend.position = "bottom",
        panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        panel.spacing.x = unit(1.5, "lines"))
dev.off()

prop <- dat_full %>% dplyr::select(cntry, thermo.diff, rur)
prop$pos.diff <- ifelse(prop$thermo.diff>0, 1, 0)

prop %>%
  count(pos.diff) %>%
  mutate(percent = n / sum(n) * 100)

prop %>%
  group_by(cntry) %>%
  count(pos.diff) %>%
  mutate(percent = n / sum(n) * 100)

prop %>%
  group_by(rur) %>%
  count(pos.diff) %>%
  mutate(percent = n / sum(n) * 100)

grp.mean_wide <- spread(grp.mean, rur, Mean)
grp.mean_wide$diff <- grp.mean_wide[,3]-grp.mean_wide[,2]

# 5.2 Figure 2 #############################################################
thermo1 <- lm(thermo.diff ~ urban_move + cntry, data = subset(dat_full, dat_full$rur=="Urban residence"))
summary(thermo1)
dat_thermo1 <- broom::tidy(thermo1)
dat_thermo1 %<>% filter(term=="urban_moveYes")
dat_thermo1$sample <- "Urban residence"

thermo2 <- lm(thermo.diff ~ rural_move + cntry, data = subset(dat_full, dat_full$rur=="Rural residence"))
summary(thermo2)
dat_thermo2 <- broom::tidy(thermo2)
dat_thermo2 %<>% filter(term=="rural_moveYes")
dat_thermo2$sample <- "Rural residence"

dat_thermo_fin <- rbind(dat_thermo1, dat_thermo2)

dat_thermo_fin$term[dat_thermo_fin$term=="urban_moveYes"] <- "Urban move \n(Rural → Urban)"
dat_thermo_fin$term[dat_thermo_fin$term=="rural_moveYes"] <- "Rural move \n(Urban → Rural)"

dat_thermo_fin$sample <- as.factor(dat_thermo_fin$sample)
dat_thermo_fin$sample <- fct_rev(dat_thermo_fin$sample)
dat_thermo_fin$estimate <- as.numeric(dat_thermo_fin$estimate)
dat_thermo_fin$std.error <- as.numeric(dat_thermo_fin$std.error)
dat_thermo_fin$var <- dat_thermo_fin$term
dat_thermo_fin$var <- factor(as.factor(dat_thermo_fin$var), levels = c("Urban move \n(Rural → Urban)", "Rural move \n(Urban → Rural)"))

Fig2 <- ggplot(dat_thermo_fin, aes(x = var, y = estimate, ))
Fig2 <- Fig2 + geom_hline(yintercept = 0, color = "black", lty = 2)

Fig2 <- Fig2 + geom_linerange(aes(x = var, ymin = estimate - std.error*interval1,
                                      ymax = estimate + std.error*interval1, color = sample),
                                  lwd = 1, position = position_dodge(width = 1/2))  +
  scale_color_manual(name = "Subsample:", values = c("black", "red")) 
  
Fig2 <- Fig2 + geom_pointrange(aes(x = var, y = estimate, ymin = estimate - std.error*interval2,
                                       ymax = estimate + std.error*interval2, color = sample, shape = sample),
                                   lwd = 1/2, position = position_dodge(width = 1/2)) +
  scale_color_manual(name = "Subsample:", values = c("black", "red")) +
  scale_shape_manual(name = "Subsample:", values = c(1, 2))

Fig2 <- Fig2 + xlab("Independent variable:") + ylab("Estimated coefficient \n Dependent variable: place-based affective polarisation") +  
  theme_pubr() +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_text(), panel.spacing.x = unit(1, "lines")) 

Fig2 <- Fig2 + coord_flip(ylim = c(-4, 4)) 

png(file=here("Figure2.png"),
    width = 8, height = 3, units = 'in', res = 800)
Fig2 
dev.off()

# 5.3 Figure 3 #############################################################
thermo1 <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence"))
summary(thermo1)
dat_thermo1 <- broom::tidy(thermo1)
dat_thermo1 %<>% filter(term=="resent_munis.std" | term=="id.std")
dat_thermo1$sample <- "Urban residence"

thermo2 <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence"))
summary(thermo2)
dat_thermo2 <- broom::tidy(thermo2)
dat_thermo2 %<>% filter(term=="resent_munis.std" | term=="id.std")
dat_thermo2$sample <- "Rural residence"

dat_thermo_fin <- rbind(dat_thermo1, dat_thermo2)

dat_thermo_fin$term[dat_thermo_fin$term=="resent_munis.std"] <- "Place-based \nresentment (Std.)"
dat_thermo_fin$term[dat_thermo_fin$term=="id.std"] <- "Place-based \nidentity (Std.)"

dat_thermo_fin$sample <- as.factor(dat_thermo_fin$sample)
dat_thermo_fin$sample <- fct_rev(dat_thermo_fin$sample)
dat_thermo_fin$estimate <- as.numeric(dat_thermo_fin$estimate)
dat_thermo_fin$std.error <- as.numeric(dat_thermo_fin$std.error)
dat_thermo_fin$var <- dat_thermo_fin$term
dat_thermo_fin$var <- factor(as.factor(dat_thermo_fin$var), levels = c("Place-based \nidentity (Std.)", "Place-based \nresentment (Std.)"))


Fig3 <- ggplot(dat_thermo_fin, aes(x = var, y = estimate, ))
Fig3 <- Fig3 + geom_hline(yintercept = 0, color = "black", lty = 2)

Fig3 <- Fig3 + geom_linerange(aes(x = var, ymin = estimate - std.error*interval1,
                                    ymax = estimate + std.error*interval1, color = sample),
                                lwd = 1, position = position_dodge(width = 1/2))  +
  scale_color_manual(name = "Subsample:", values = c("black", "red"))

Fig3 <- Fig3 + geom_pointrange(aes(x = var, y = estimate, ymin = estimate - std.error*interval2,
                                     ymax = estimate + std.error*interval2, color = sample, shape = sample),
                                 lwd = 1/2, position = position_dodge(width = 1/2)) +
  scale_color_manual(name = "Subsample:", values = c("black", "red")) +
  scale_shape_manual(name = "Subsample:", values = c(1, 2))

Fig3 <- Fig3 + xlab("Independent variable:") + ylab("Estimated coefficient \n Dependent variable: place-based affective polarisation") +  
  theme_pubr() +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_text(), panel.spacing.x = unit(1, "lines")) 

Fig3 <- Fig3 + coord_flip(ylim = c(-3.5, 10.5)) 


png(file=here("Figure3.png"),
    width = 8, height = 3, units = 'in', res = 800)
Fig3 
dev.off()

# 5.4 Figure 4 #############################################################
thermo_in1 <- lm(thermo.in ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence"))
dat_in_urb <- broom::tidy(thermo_in1)
dat_in_urb %<>% filter(term=="resent_munis.std" | term=="id.std")
dat_in_urb$sample1 <- "Subsample:\n Urban residence"
dat_in_urb$sample2 <- "In-group affect"

thermo_out1 <- lm(thermo.out ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence"))
dat_out_urb <- broom::tidy(thermo_out1)
dat_out_urb %<>% filter(term=="resent_munis.std" | term=="id.std")
dat_out_urb$sample1 <- "Subsample:\n Urban residence"
dat_out_urb$sample2 <- "Out-group affect"

thermo_in2 <- lm(thermo.in ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence"))
dat_in_rur <- broom::tidy(thermo_in2)
dat_in_rur %<>% filter(term=="resent_munis.std" | term=="id.std")
dat_in_rur$sample1 <- "Subsample:\n Rural residence"
dat_in_rur$sample2 <- "In-group affect"

thermo_out2 <- lm(thermo.out ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence"))
dat_out_rur <- broom::tidy(thermo_out2)
dat_out_rur %<>% filter(term=="resent_munis.std" | term=="id.std")
dat_out_rur$sample1 <- "Subsample:\n Rural residence"
dat_out_rur$sample2 <- "Out-group affect"

dat_inout_fin1 <- rbind(dat_in_urb, dat_out_urb)
dat_inout_fin2 <- rbind(dat_inout_fin1, dat_in_rur)
dat_inout_fin3 <- rbind(dat_inout_fin2, dat_out_rur)
dat_inout_fin <- dat_inout_fin3

dat_inout_fin$term[dat_inout_fin$term=="resent_munis.std"] <- "Place-based \nresentment (Std.)"
dat_inout_fin$term[dat_inout_fin$term=="id.std"] <- "Place-based \nidentity (Std.)"


dat_inout_fin$term <- as.factor(dat_inout_fin$term)
dat_inout_fin$sample1 <- as.factor(dat_inout_fin$sample1)
dat_inout_fin$sample1 <- factor(dat_inout_fin$sample1, levels = c("Subsample:\n Urban residence", "Subsample:\n Rural residence"))
dat_inout_fin$sample2 <- as.factor(dat_inout_fin$sample2)
dat_inout_fin$sample2 <- fct_rev(dat_inout_fin$sample2)
dat_inout_fin$estimate <- as.numeric(dat_inout_fin$estimate)
dat_inout_fin$std.error <- as.numeric(dat_inout_fin$std.error)

Fig4 <- ggplot(dat_inout_fin, aes(x = term, y = estimate))
Fig4 <- Fig4 + geom_hline(yintercept = 0, color = "black", lty = 2)

Fig4 <- Fig4 + geom_linerange(aes(x = term, ymin = estimate - std.error*interval1,
                                  ymax = estimate + std.error*interval1, color = sample2),
                              lwd = 1, position = position_identity()) +
  scale_color_manual(values = c("black", "red"), labels=c("Out-group affect", "In-group affect"),
                     name = "Dependent variable:") + 
  scale_shape_manual(values = c(1, 2), labels=c("Out-group affect", "In-group affect"),
                     name = "Dependent variable:") 

Fig4 <- Fig4 + geom_pointrange(aes(x = term, y = estimate, ymin = estimate - std.error*interval2,
                                   ymax = estimate + std.error*interval2, color = sample2, shape = sample2),
                               lwd = 1/2, position = position_identity()) +
  scale_color_manual(values = c("black", "red"), labels=c("Out-group affect", "In-group affect"),
                     name = "Dependent variable:") + 
  scale_shape_manual(values = c(1, 2), labels=c("Out-group affect", "In-group affect"),
                     name = "Dependent variable:") 

Fig4 <- Fig4 + xlab("Independent variable:") + ylab("Estimated coefficient") + 
  theme_pubr() +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        panel.spacing.x = unit(1, "lines")) 

png(file=here("Figure4.png"),
    width = 8, height = 2.75, units = 'in', res = 800)
Fig4 + facet_wrap(vars(sample1)) + coord_flip(ylim = c(-5, 10)) 
dev.off()

# 5.5 Figure 5 #############################################################
GALTAN <- lm(galtan_vote ~ thermo.diff*rur + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = dat_full)

GALTAN <- predict_response(GALTAN, terms = c("thermo.diff", "rur"))

Fig5_1 <- ggplot(GALTAN, aes(x = x, y = predicted, group = group)) +
  geom_line(aes(linetype=group)) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill=group), alpha = .2) +
  labs(y = "← GAL-TAN →",
       x = "Thermometer differential")  +
  scale_fill_manual(values = c("black", "grey")) + 
  theme_pubr() +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_blank(), panel.spacing.x = unit(1.5, "lines"))

hist_theme <- theme_void() + theme(legend.position = "none")
dat_hist <- dat_full %>% dplyr::select(galtan_vote, thermo.diff, rur, gender, age.std, educ, inc_dec, migbck, lr.std, cntry) %>% na.omit()

Fig5_2 <- ggplot(dat_hist, aes(x = thermo.diff, fill = rur, linetype=rur)) + 
  geom_density(alpha=.5) +
  scale_fill_manual(values = c("black", "grey")) +  
  hist_theme

png(file=here("Figure5.png"),
    width = 8, height = 6, units = 'in', res = 800)
Fig5_2 + 
  Fig5_1 + plot_layout(byrow = FALSE, 
                        widths = c(1, 1), 
                        heights = unit(c(3, 1), c('cm', 'null')), ncol = 1)
dev.off()

# 5.6 Figure 6 #############################################################
png(file=here("Figure6.png"),
    width = 8, height = 6, units = 'in', res = 800)
ggscatter(CHES_ts, x="urban_rural", y="galtan",
          add = "reg.line",
          add.params = list(color = "blue", fill = "lightgray"), 
          conf.int = TRUE) +
  labs(x = "← Urban-Rural →", y = "← GAL-TAN →") + 
  scale_x_continuous(breaks = c(0.0, 2.5, 5, 7.5, 10.0), limits = c(-0.5,10.5)) +
  scale_y_continuous(breaks = c(0.0, 2.5, 5, 7.5, 10.0), limits = c(-3.75,13.75)) +
  geom_text_repel(aes(label = party), color = "black", size = 3.5, max.overlaps = Inf) +
  facet_wrap(CHES_ts$country, ncol =3) +
  stat_cor(method = "pearson") +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_blank(), panel.spacing.x = unit(1.5, "lines")) + 
  coord_axes_inside(xintercept = 5, yintercept = 5) 
dev.off()

# 5.7 Figure 7 #############################################################
GALTAN <- lm(galtan_vote ~ thermo.in*rur + thermo.out*rur + gender + age.std + educ + inc_dec + migbck + lr.std + cntry,
            data = dat_full)
summary(GALTAN)

GALTAN <- predict_response(GALTAN, terms = c("thermo.in", "rur"))


Fig7_1 <- ggplot(GALTAN, aes(x = x, y = predicted, group = group)) +
  geom_line(aes(linetype=group)) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill=group), alpha = .2) +
  labs(y = "← GAL-TAN →",
       x = "In-group affect")  +
  ylim(4.5, 8) + 
  scale_fill_manual(values = c("black", "grey")) + 
  theme_pubr() +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_blank(), panel.spacing.x = unit(1.5, "lines"))

hist_theme <- theme_void() + theme(legend.position = "none")
dat_hist <- dat_full %>% dplyr::select(galtan_vote, thermo.in, rur, gender, age.std, educ, inc_dec, migbck, lr.std, cntry) %>% na.omit()

Fig7_2 <- ggplot(dat_hist, aes(x = thermo.in, fill = rur, linetype=rur)) + 
  geom_density(alpha=.5) +
  scale_fill_manual(values = c("black", "grey")) +  
  hist_theme

GALTAN <- lm(galtan_vote ~ thermo.in*rur + thermo.out*rur + gender + age.std + educ + inc_dec + migbck + lr.std + cntry,
             data = dat_full)

GALTAN <- predict_response(GALTAN, terms = c("thermo.out", "rur"))

Fig7_3 <- ggplot(GALTAN, aes(x = x, y = predicted, group = group)) +
  geom_line(aes(linetype=group)) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill=group), alpha = .2) +
  labs(y = "← GAL-TAN →",
       x = "Out-group affect")  +
  ylim(4.5, 8) + 
  scale_fill_manual(values = c("black", "grey")) + 
  theme_pubr() +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_blank(), panel.spacing.x = unit(1.5, "lines"))
Fig7_3

hist_theme <- theme_void() + theme(legend.position = "none")
dat_hist <- dat_full %>% dplyr::select(galtan_vote, thermo.out, rur, gender, age.std, educ, inc_dec, migbck, lr.std, cntry) %>% na.omit()

Fig7_4 <- ggplot(dat_hist, aes(x = thermo.out, fill = rur, linetype=rur)) + 
  geom_density(alpha=.5) +
  scale_fill_manual(values = c("black", "grey")) +  
  hist_theme

png(file=here("Figure7.png"),
    width = 8, height = 6, units = 'in', res = 800)
Fig7_2 + Fig7_1 + Fig7_4 + Fig7_3 + plot_layout(byrow = FALSE, 
                                        widths = c(1, 1), 
                                        heights = unit(c(3, 1), c('cm', 'null')))
dev.off()

# 6. Appendix ##############################################################
# 6.1 Table A1 #############################################################
dat_summary <- dat_full %>% dplyr::select(gender, age, age.std, educ, inc_dec, 
                                   galtan_vote, radicalright, green,
                                   thermo.diff, thermo.diff.std,
                                   thermo.in, thermo.in.std,
                                   thermo.out, thermo.out.std,
                                   trait.diff, trait.diff.std,
                                   trait.in, trait.out,
                                   int.in, open.in, hon.in, self.in, hypo.in,
                                   int.out, open.out, hon.out, self.out, hypo.out,
                                   resent_munis, resent_munis.std,
                                   eco, rep1, rep2, cult1, cult2,
                                   id, id.std,
                                   urid_urb,
                                   rurid_rur,
                                   lr, lr.std,
                                   plid, plid.std,
                                   rur, rur_cont,
                                   migbck,
                                   immi, immi.std,
                                   res_num, resbio_num, 
                                   movement, urban_move, rural_move,
                                   reg_unemp, pop_delta, popdens_n3, EQI,
                                   reg_unemp.std, pop_delta.std)

dat_summary$gender.num <- as.numeric(dat_summary$gender)
dat_summary$gender.num <- ifelse(dat_summary$gender.num==1, 0, 1)

dat_summary$educ.num <- as.numeric(dat_summary$educ)
dat_summary$educ.num <- ifelse(dat_summary$educ.num==1, 0, 1)

dat_summary$rur.num <- as.numeric(dat_summary$rur)
dat_summary$rur.num <- ifelse(dat_summary$rur.num==1, 0, 1)

dat_summary$rur_cont <- as.character(dat_summary$rur_cont)
dat_summary$rur_rec <- recode(dat_summary$rur_cont, "Very urban"="0", "Rather urban"="1",
                              "Rather rural"="2", "Very rural"="3")
dat_summary$rur_cont <- as.numeric(dat_summary$rur_rec)

dat_summary$migbck.num <- as.numeric(dat_summary$migbck)
dat_summary$migbck.num <- ifelse(dat_summary$migbck.num==1, 0, 1)

dat_summary$urban_move.num <- as.numeric(dat_summary$urban_move)
dat_summary$urban_move.num <- ifelse(dat_summary$urban_move.num==1, 0, 1)

dat_summary$rural_move.num <- as.numeric(dat_summary$rural_move)
dat_summary$rural_move.num <- ifelse(dat_summary$rural_move.num==1, 0, 1)

stargazer(dat_summary, digits=2, type = "text")

# 6.2 Figure A1 ############################################################
grp.mean <- dat_full %>%
  group_by(cntry, rur) %>%
  dplyr::summarize(Mean = mean(trait.diff, na.rm=TRUE))

png(file=here("FigureA1.png"),
    width = 8, height = 6, units = 'in', res = 800)
ggplot(dat_full, aes(x = trait.diff, fill = rur)) + 
  geom_density(alpha=.5) +
  scale_fill_manual(name = "Subsample:", values = c("black", "grey")) + 
  labs(y = "Density",
       x = "Trait-rating differential") + 
  facet_wrap(vars(cntry)) +
  geom_vline(data=grp.mean, aes(xintercept=Mean, color = rur),
             linetype="dashed") + 
  scale_color_manual(name = "Subsample:", values = c("black", "grey")) + 
  theme_pubr() +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_text(), panel.spacing.x = unit(1.5, "lines"))
dev.off()

# 6.3 Figure A2 ############################################################
thermo1 <- lm(trait.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence"))
summary(thermo1)
dat_thermo1 <- broom::tidy(thermo1)
dat_thermo1 %<>% filter(term=="resent_munis.std" | term=="id.std")
dat_thermo1$sample <- "Urban residence"


thermo2 <- lm(trait.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence"))
summary(thermo2)
dat_thermo2 <- broom::tidy(thermo2)
dat_thermo2 %<>% filter(term=="resent_munis.std" | term=="id.std")
dat_thermo2$sample <- "Rural residence"

dat_thermo_fin <- rbind(dat_thermo1, dat_thermo2)

dat_thermo_fin$term[dat_thermo_fin$term=="resent_munis.std"] <- "Place-based \nresentment (Std.)"
dat_thermo_fin$term[dat_thermo_fin$term=="id.std"] <- "Place-based \nidentity (Std.)"

dat_thermo_fin$sample <- as.factor(dat_thermo_fin$sample)
dat_thermo_fin$sample <- fct_rev(dat_thermo_fin$sample)
dat_thermo_fin$estimate <- as.numeric(dat_thermo_fin$estimate)
dat_thermo_fin$std.error <- as.numeric(dat_thermo_fin$std.error)
dat_thermo_fin$var <- dat_thermo_fin$term
dat_thermo_fin$var <- factor(as.factor(dat_thermo_fin$var), levels = c("Place-based \nidentity (Std.)", "Place-based \nresentment (Std.)"))

FigA2 <- ggplot(dat_thermo_fin, aes(x = var, y = estimate, ))
FigA2 <- FigA2 + geom_hline(yintercept = 0, color = "black", lty = 2)

FigA2 <- FigA2 + geom_linerange(aes(x = var, ymin = estimate - std.error*interval1,
                                      ymax = estimate + std.error*interval1, color = sample),
                                  lwd = 1, position = position_dodge(width = 1/2))  +
  scale_color_manual(name = "Subsample:", values = c("black", "red")) 

FigA2 <- FigA2 + geom_pointrange(aes(x = var, y = estimate, ymin = estimate - std.error*interval2,
                                       ymax = estimate + std.error*interval2, color = sample, shape = sample),
                                   lwd = 1/2, position = position_dodge(width = 1/2)) +
  scale_color_manual(name = "Subsample:", values = c("black", "red")) +
  scale_shape_manual(name = "Subsample:", values = c(1, 2))

FigA2 <- FigA2 + xlab("Independent variable:") + ylab("Estimated coefficient \n Dependent variable: place-based affective polarisation") +  
  theme_pubr() +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_text(), panel.spacing.x = unit(1, "lines")) 

FigA2 <- FigA2 + coord_flip(ylim = c(0, 0.5)) 

png(file=here("FigureA2.png"),
    width = 8, height = 3, units = 'in', res = 800)
FigA2 
dev.off()

# 6.4 Table A5 #############################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), ncol = 7)) 
attr(FE_add, "position") = 19

TableA5 <- list("(1)" <- lm(trait.diff ~ resent_munis.std + id.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                "(2)" <- lm(trait.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                "(3)" <- lm(trait.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                
                "(4)" <- lm(trait.diff ~ resent_munis.std + id.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence")),
                "(5)" <- lm(trait.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + cntry, data = subset(dat_full, dat_full$rur=="Rural residence")),
                "(6)" <- lm(trait.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence"))
)

modelsummary(TableA5, fmt = 3, stars = T,
             title = 'OLS regression results: place-based affective polarisation on place-based resentment and place-based identity, by self-classified urban-rural residence (trait-rating differential).', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("resent_munis.std" = "Place-based resentment (Std.)",
                          "id.std" = "Place-based identity (Std.)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "migbckYes" = "Migration background (b.=no)",
                          "lr.std" = "Left-right (Std.)",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.5 Figure A3 ############################################################
GALTAN <- lm(galtan_vote ~ trait.diff*rur + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = dat_full)

GALTAN <- predict_response(GALTAN, terms = c("trait.diff", "rur"))

FigA3.1 <- ggplot(GALTAN, aes(x = x, y = predicted, group = group)) +
  geom_line(aes(linetype=group)) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill=group), alpha = .2) +
  labs(y = "← GAL-TAN →",
       x = "Trait-rating differential")  +
  scale_fill_manual(values = c("black", "grey")) + 
  theme_pubr() +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_blank(), panel.spacing.x = unit(1.5, "lines"))
FigA3.1

hist_theme <- theme_void() + theme(legend.position = "none")
dat_hist <- dat_full %>% dplyr::select(galtan_vote, trait.diff, rur, gender, age.std, educ, inc_dec, migbck, lr.std, cntry) %>% na.omit()

FigA3.2 <- ggplot(dat_hist, aes(x = trait.diff, fill = rur, linetype=rur)) + 
  geom_density(alpha=.5) +
  scale_fill_manual(values = c("black", "grey")) +  
  hist_theme

png(file=here("FigureA3.png"),
    width = 8, height = 6, units = 'in', res = 800)
FigA3.2 + 
  FigA3.1 + plot_layout(byrow = FALSE, 
                     widths = c(1, 1), 
                     heights = unit(c(3, 1), c('cm', 'null')), ncol = 1)
dev.off()

# 6.6 Table A6 #############################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes"), ncol = 4)) 
attr(FE_add, "position") = 21

TableB6 <- list("(1)" <- lm(galtan_vote ~ trait.diff.std*rur + cntry, data = dat_full),
                "(2)" <- lm(galtan_vote ~ trait.diff.std*rur + gender + age.std + educ + inc_dec + migbck + cntry, data = dat_full),
                "(3)" <- lm(galtan_vote ~ trait.diff.std*rur + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = dat_full)
)

modelsummary(TableB6, fmt = 3, stars = T,
             title = 'OLS regression results: GAL-TAN voting on place-based resentment (trait-rating differential) .', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("trait.diff.std" = "Trait-rating differential (Std.)",
                          "rurRural residence" = "Rural residence (b.=urban residence)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "migbckYes" = "Migration background (b.=no)",
                          "lr.std" = "Left-right (Std.)",
                          "trait.diff.std:rurRural residence" = "Trait-rating differential (Std.) X Rural residence",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.7 Figure A4 ############################################################
grp.mean <- dat_full %>%
  group_by(cntry, rur) %>%
  dplyr::summarize(Mean = mean(resent_munis, na.rm=TRUE))

png(file=here("FigureA4.png"),
    width = 8, height = 6, units = 'in', res = 800)
ggplot(dat_full, aes(x = resent_munis, fill = rur)) + 
  geom_density(alpha=.5) +
  scale_fill_manual(name = "Subsample:", values = c("black", "grey")) + 
  labs(y = "Density",
       x = "Place-based resentment") + 
  facet_wrap(vars(cntry)) +
  geom_vline(data=grp.mean, aes(xintercept=Mean, color = rur),
             linetype="dashed") + 
  scale_color_manual(name = "Subsample:", values = c("black", "grey")) + 
  theme_pubr() +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_text(), panel.spacing.x = unit(1.5, "lines"))
dev.off()

# 6.8 Figure A5 ############################################################
resent.cor.dat <- dat_full[,c("rur","resent_munis" ,"eco", "rep1","rep2","cult1","cult2","thermo.diff")]
names(resent.cor.dat) <- c("rur","Place-based resentment", 
                           "Economic", "Representation (A)", "Representation (B)", 
                           "Culture (A)","Culture (B)" ,"Thermometer differential")

resent.cor.dat.urban <- subset(resent.cor.dat, resent.cor.dat$rur=="Urban residence")
resent.cor.dat.urban$rur <- NULL 
resent.cor.dat.rural <- subset(resent.cor.dat, resent.cor.dat$rur=="Rural residence")
resent.cor.dat.rural$rur <- NULL 

corr.urban <- round(cor(resent.cor.dat.urban, use = "complete.obs"), 2)
corr.urban.df <- as.data.frame(as.table(corr.urban))
corr.urban.df$rur <- "Subsample:\n Urban residence"

corr.rural <- round(cor(resent.cor.dat.rural, use = "complete.obs"), 2)
corr.rural.df <- as.data.frame(as.table(corr.rural))
corr.rural.df$rur <- "Subsample:\n Rural residence"

corr.df <- rbind(corr.urban.df, corr.rural.df)
corr.df$rur <- as.factor(corr.df$rur)
corr.df <- within(corr.df, rur <- relevel(rur, ref = "Subsample:\n Urban residence"))

png(file=here("FigureA5.png"),
    width = 8, height = 6, units = 'in', res = 800)
ggplot(corr.df, aes(Var1, Var2, fill = Freq)) +
  geom_tile(color = "white") +
  geom_text(aes(label = round(Freq, 2)), color = "black", size = 3) +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, 
                       limit = c(-1, 1), space = "Lab", name = "Correlation") +
  theme_minimal() +
  facet_wrap(~rur) +  
  labs(
    x = "", y = "") +
  theme_pubr() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),       
        legend.title=element_blank(), panel.spacing.x = unit(1.5, "lines"))
dev.off()

# 6.9 Figure A6 ############################################################
grp.mean <- dat_full %>%
  group_by(cntry, rur) %>%
  dplyr::summarize(Mean = mean(id, na.rm=TRUE))

png(file=here("FigureA6.png"),
    width = 8, height = 6, units = 'in', res = 800)
ggplot(dat_full, aes(x = id, fill = rur)) + 
  geom_density(alpha=.5) +
  scale_fill_manual(name = "Subsample:", values = c("black", "grey")) + 
  labs(y = "Density",
       x = "Place-based identity") + 
  facet_wrap(vars(cntry)) +
  geom_vline(data=grp.mean, aes(xintercept=Mean, color = rur),
             linetype="dashed") + 
  scale_color_manual(name = "Subsample:", values = c("black", "grey")) + 
  theme_pubr() +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_text(), panel.spacing.x = unit(1.5, "lines"))
dev.off()

# 6.10 Figure A7 ###########################################################
grp.mean <- dat_full %>%
  group_by(cntry, rur) %>%
  dplyr::summarize(Mean = mean(plid, na.rm=TRUE))

png(file=here("FigureA7.png"),
    width = 8, height = 6, units = 'in', res = 800)
ggplot(dat_full, aes(x = plid, fill = rur)) + 
  geom_density(alpha=.5) +
  scale_fill_manual(name = "Subsample:", values = c("black", "grey")) + 
  labs(y = "Density",
       x = "Attachment to place of residence") + 
  facet_wrap(vars(cntry)) +
  geom_vline(data=grp.mean, aes(xintercept=Mean, color = rur),
             linetype="dashed") + 
  scale_color_manual(name = "Subsample:", values = c("black", "grey")) + 
  theme_pubr() +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_text(), panel.spacing.x = unit(1.5, "lines"))
dev.off()

# 6.11 Figure A8 ############################################################
thermo1 <- lm(thermo.diff ~ resent_munis.std + plid.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence"))
summary(thermo1)
dat_thermo1 <- broom::tidy(thermo1)
dat_thermo1 %<>% filter(term=="resent_munis.std" | term=="plid.std")
dat_thermo1$sample <- "Urban residence"


thermo2 <- lm(thermo.diff ~ resent_munis.std + plid.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence"))
summary(thermo2)
dat_thermo2 <- broom::tidy(thermo2)
dat_thermo2 %<>% filter(term=="resent_munis.std" | term=="plid.std")
dat_thermo2$sample <- "Rural residence"

dat_thermo_fin <- rbind(dat_thermo1, dat_thermo2)

dat_thermo_fin$term[dat_thermo_fin$term=="resent_munis.std"] <- "Place-based \nresentment (Std.)"
dat_thermo_fin$term[dat_thermo_fin$term=="plid.std"] <- "Attachment to \nplace of residence (Std.)"

dat_thermo_fin$sample <- as.factor(dat_thermo_fin$sample)
dat_thermo_fin$sample <- fct_rev(dat_thermo_fin$sample)
dat_thermo_fin$estimate <- as.numeric(dat_thermo_fin$estimate)
dat_thermo_fin$std.error <- as.numeric(dat_thermo_fin$std.error)
dat_thermo_fin$var <- dat_thermo_fin$term
dat_thermo_fin$var <- factor(as.factor(dat_thermo_fin$var), levels = c("Attachment to \nplace of residence (Std.)", "Place-based \nresentment (Std.)"))

FigA8 <- ggplot(dat_thermo_fin, aes(x = var, y = estimate, ))
FigA8 <- FigA8 + geom_hline(yintercept = 0, color = "black", lty = 2)

FigA8 <- FigA8 + geom_linerange(aes(x = var, ymin = estimate - std.error*interval1,
                                    ymax = estimate + std.error*interval1, color = sample),
                                lwd = 1, position = position_dodge(width = 1/2))  +
  scale_color_manual(name = "Subsample:", values = c("black", "red"))

FigA8 <- FigA8 + geom_pointrange(aes(x = var, y = estimate, ymin = estimate - std.error*interval2,
                                     ymax = estimate + std.error*interval2, color = sample, shape = sample),
                                 lwd = 1/2, position = position_dodge(width = 1/2)) +
  scale_color_manual(name = "Subsample:", values = c("black", "red")) +
  scale_shape_manual(name = "Subsample:", values = c(1, 2))

FigA8 <- FigA8 + xlab("Independent variable:") + ylab("Estimated coefficient \n Dependent variable: place-based affective polarisation") +  
  theme_pubr() +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_text(), panel.spacing.x = unit(1, "lines")) 

FigA8 <- FigA8 + coord_flip(ylim = c(0, 10.5)) 

png(file=here("FigureA8.png"),
    width = 8, height = 3, units = 'in', res = 800)
FigA8 
dev.off()

# 6.12 Table A8 ############################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), ncol = 7)) 
attr(FE_add, "position") = 19

TableA8 <- list("(1)" <- lm(thermo.diff ~ resent_munis.std + plid.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                "(2)" <- lm(thermo.diff ~ resent_munis.std + plid.std + gender + age.std + educ + inc_dec + migbck + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                "(3)" <- lm(thermo.diff ~ resent_munis.std + plid.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                
                "(4)" <- lm(thermo.diff ~ resent_munis.std + plid.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence")),
                "(5)" <- lm(thermo.diff ~ resent_munis.std + plid.std + gender + age.std + educ + inc_dec + migbck + cntry, data = subset(dat_full, dat_full$rur=="Rural residence")),
                "(6)" <- lm(thermo.diff ~ resent_munis.std + plid.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence"))
)

modelsummary(TableA8, fmt = 3, stars = T,
             title = 'OLS regression results: place-based affective polarisation on place-based resentment and attachment to place of residence, by self-classified urban-rural residence (trait-rating differential).', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("resent_munis.std" = "Place-based resentment (Std.)",
                          "plid.std" = "Attachment to place of residence (Std.)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "migbckYes" = "Migration background (b.=no)",
                          "lr.std" = "Left-right (Std.)",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.13 Figure A9 ###########################################################
png(file=here("FigureA9.png"),
    width = 8, height = 6, units = 'in', res = 800)
ggplot(dat_full, aes(x = galtan_vote)) + 
  geom_histogram(alpha=.5) +
  labs(y = "Count",
       x = "← GAL-TAN →") + 
  facet_wrap(vars(cntry)) +
  scale_color_manual(values = c("grey", "black")) +
  theme_pubr() +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_blank(), panel.spacing.x = unit(1.5, "lines"))
dev.off()

# 6.14 Figure A10 ##########################################################
GALTAN <- lm(radicalright ~ thermo.diff*rur + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = dat_full)

GALTAN <- predict_response(GALTAN, terms = c("thermo.diff", "rur"))

FigA10.1 <- ggplot(GALTAN, aes(x = x, y = predicted, group = group)) +
  geom_line(aes(linetype=group)) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill=group), alpha = .2) +
  labs(y = "Pr(Radical right vote)",
       x = "Thermometer differential")  +
  scale_fill_manual(values = c("black", "grey")) + 
  theme_pubr() +
  theme(legend.position="none", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_blank(), panel.spacing.x = unit(1.5, "lines"))
FigA10.1

hist_theme <- theme_void() + theme(legend.position = "none")
dat_hist <- dat_full %>% dplyr::select(radicalright, thermo.diff, rur, gender, age.std, educ, inc_dec, migbck, lr.std, cntry) %>% na.omit()

FigA10.2 <- ggplot(dat_hist, aes(x = thermo.diff, fill = rur, linetype=rur)) + 
  geom_density(alpha=.5) +
  scale_fill_manual(values = c("black", "grey"), guide ="none") +
  scale_linetype_manual(values=c("solid", "dashed"), guide ="none")+
  hist_theme

GALTAN <- lm(green ~ thermo.diff*rur + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = dat_full)

GALTAN <- predict_response(GALTAN, terms = c("thermo.diff", "rur"))

FigA10.3 <- ggplot(GALTAN, aes(x = x, y = predicted, group = group)) +
  geom_line(aes(linetype=group)) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill=group), alpha = .2) +
  labs(y = "Pr(Green vote)",
       x = "Thermometer differential")  +
  scale_fill_manual(values = c("black", "grey")) + 
  theme_pubr() +
  theme(legend.position="none", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_blank(), panel.spacing.x = unit(1.5, "lines"))

hist_theme <- theme_void() + theme(legend.position = "none")
dat_hist <- dat_full %>% dplyr::select(green, thermo.diff, rur, gender, age.std, educ, inc_dec, migbck, lr.std, cntry) %>% na.omit()

FigA10.4 <- ggplot(dat_hist, aes(x = thermo.diff, fill = rur, linetype=rur)) + 
  geom_density(alpha=.5) +
  scale_fill_manual(values = c("black", "grey"), guide ="none") +
  scale_linetype_manual(values=c("solid", "dashed"), guide ="none")+
  hist_theme

FigA10 <- FigA10.2 + FigA10.1 + FigA10.4 + FigA10.3 + plot_layout(byrow = FALSE, 
                                                                  widths = c(1, 1), 
                                                                  heights = unit(c(3, 1), c('cm', 'null')))

png(file=here("FigureA10.png"),
    width = 8, height = 6, units = 'in', res = 800)
FigA10 + plot_layout(guides = "collect") & theme(legend.position = 'bottom')
dev.off()

# 6.15 Table A9 #############################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), ncol = 7)) 
attr(FE_add, "position") = 21

TableA9 <- list("(1)" <- lm(radicalright ~ thermo.diff.std*rur + cntry, data = dat_full),
                "(2)" <- lm(radicalright ~ thermo.diff.std*rur + gender + age.std + educ + inc_dec + migbck + cntry, data = dat_full),
                "(3)" <- lm(radicalright ~ thermo.diff.std*rur + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = dat_full),
                
                "(4)" <- lm(green ~ thermo.diff.std*rur + cntry, data = dat_full),
                "(5)" <- lm(green ~ thermo.diff.std*rur + gender + age.std + educ + inc_dec + migbck + cntry, data = dat_full),
                "(6)" <- lm(green ~ thermo.diff.std*rur + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = dat_full)
)

modelsummary(TableA9, fmt = 3, stars = T,
             title = 'OLS regression results: radical right and green voting on place-based resentment.', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("thermo.diff.std" = "Thermometer differential (Std.)",
                          "rurRural residence" = "Rural residence (b.=urban residence)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "migbckYes" = "Migration background (b.=no)",
                          "lr.std" = "Left-right (Std.)",
                          "thermo.diff.std:rurRural residence" = "Thermometer differential (Std.) X Rural residence",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.16 Table A10 ############################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), ncol = 7)) 
attr(FE_add, "position") = 21

TableA10 <- list("(1)" <- glm(radicalright ~ thermo.diff.std*rur + cntry, data = dat_full, family = "binomial"),
                "(2)" <- glm(radicalright ~ thermo.diff.std*rur + gender + age.std + educ + inc_dec + migbck + cntry, data = dat_full, family = "binomial"),
                "(3)" <- glm(radicalright ~ thermo.diff.std*rur + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = dat_full, family = "binomial"),
                
                "(4)" <- glm(green ~ thermo.diff.std*rur + cntry, data = dat_full, family = "binomial"),
                "(5)" <- glm(green ~ thermo.diff.std*rur + gender + age.std + educ + inc_dec + migbck + cntry, data = dat_full, family = "binomial"),
                "(6)" <- glm(green ~ thermo.diff.std*rur + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = dat_full, family = "binomial")
)

modelsummary(TableA10, fmt = 3, stars = T,
             title = 'Logistic regression results: radical right and green voting on place-based resentment.', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("thermo.diff.std" = "Thermometer differential (Std.)",
                          "rurRural residence" = "Rural residence (b.=urban residence)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "migbckYes" = "Migration background (b.=no)",
                          "lr.std" = "Left-right (Std.)",
                          "thermo.diff.std:rurRural residence" = "Thermometer differential (Std.) X Rural residence",
                          "(Intercept)" = "Constant"))

# 6.17 Figure A11 ##########################################################
thermo1 <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur_cont=="Very urban"))
summary(thermo1)
dat_thermo1 <- broom::tidy(thermo1)
dat_thermo1 %<>% filter(term=="resent_munis.std" | term=="id.std")
dat_thermo1$sample <- "Very urban"

thermo2 <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur_cont=="Rather urban"))
summary(thermo2)
dat_thermo2 <- broom::tidy(thermo2)
dat_thermo2 %<>% filter(term=="resent_munis.std" | term=="id.std")
dat_thermo2$sample <- "Rather urban"

thermo3 <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck+ lr.std + cntry, data = subset(dat_full, dat_full$rur_cont=="Rather rural"))
summary(thermo3)
dat_thermo3 <- broom::tidy(thermo3)
dat_thermo3 %<>% filter(term=="resent_munis.std" | term=="id.std")
dat_thermo3$sample <- "Rather rural"

thermo4 <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur_cont=="Very rural"))
summary(thermo4)
dat_thermo4 <- broom::tidy(thermo4)
dat_thermo4 %<>% filter(term=="resent_munis.std" | term=="id.std")
dat_thermo4$sample <- "Very rural"

dat_thermo_fin1 <- rbind(dat_thermo1, dat_thermo2)
dat_thermo_fin2 <- rbind(dat_thermo_fin1, dat_thermo3)
dat_thermo_fin3 <- rbind(dat_thermo_fin2, dat_thermo4)
dat_thermo_fin <- dat_thermo_fin3

dat_thermo_fin$term[dat_thermo_fin$term=="resent_munis.std"] <- "Place-based \nresentment (Std.)"
dat_thermo_fin$term[dat_thermo_fin$term=="id.std"] <- "Place-based \nidentity (Std.)"

dat_thermo_fin$sample <- as.factor(dat_thermo_fin$sample)
dat_thermo_fin$sample <- ordered(dat_thermo_fin$sample, levels = c("Very urban", "Rather urban", "Rather rural", "Very rural"))
dat_thermo_fin$estimate <- as.numeric(dat_thermo_fin$estimate)
dat_thermo_fin$std.error <- as.numeric(dat_thermo_fin$std.error)
dat_thermo_fin$var <- dat_thermo_fin$term
dat_thermo_fin$var <- factor(as.factor(dat_thermo_fin$var), levels = c("Place-based \nidentity (Std.)", "Place-based \nresentment (Std.)"))

FigA11 <- ggplot(dat_thermo_fin, aes(x = var, y = estimate, ))
FigA11 <- FigA11 + geom_hline(yintercept = 0, color = "black", lty = 2)

FigA11 <- FigA11 + geom_linerange(aes(x = var, ymin = estimate - std.error*interval1,
                                      ymax = estimate + std.error*interval1, color = sample),
                                  lwd = 1, position = position_dodge(width = 1/2))  +
  scale_color_manual(name = "Subsample:", values = c("black", "black", "red", "red"))


FigA11 <- FigA11 + geom_pointrange(aes(x = var, y = estimate, ymin = estimate - std.error*interval2,
                                       ymax = estimate + std.error*interval2, color = sample, shape = sample),
                                   lwd = 1/2, position = position_dodge(width = 1/2)) +
  scale_color_manual(name = "Subsample:", values = c("black", "black", "red", "red")) +
  scale_shape_manual(name = "Subsample:", values=c(1,4,3,2))

FigA11 <- FigA11 + xlab("Independent variable:") + ylab("Estimated coefficient \n Dependent variable: place-based affective polarisation") +  
  theme_pubr() +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_text(), panel.spacing.x = unit(1, "lines")) 

FigA11 <- FigA11 + coord_flip(ylim = c(-3.5, 10.5)) 

png(file=here("FigureA11.png"),
    width = 8, height = 3, units = 'in', res = 800)
FigA11 
dev.off()

# 6.18 Table A11 ############################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes",
                             "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), ncol = 13)) 
attr(FE_add, "position") = 19

TableA11 <- list("(1)" <- lm(thermo.diff ~ resent_munis.std + id.std + cntry, data = subset(dat_full, dat_full$rur_cont=="Very urban")),
                 "(2)" <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + cntry, data = subset(dat_full, dat_full$rur_cont=="Very urban")),
                 "(3)" <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur_cont=="Very urban")),
                 
                 "(4)" <- lm(thermo.diff ~ resent_munis.std + id.std + cntry, data = subset(dat_full, dat_full$rur_cont=="Rather urban")),
                 "(5)" <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + cntry, data = subset(dat_full, dat_full$rur_cont=="Rather urban")),
                 "(6)" <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur_cont=="Rather urban")),
                 
                 "(7)" <- lm(thermo.diff ~ resent_munis.std + id.std + cntry, data = subset(dat_full, dat_full$rur_cont=="Rather rural")),
                 "(8)" <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + cntry, data = subset(dat_full, dat_full$rur_cont=="Rather rural")),
                 "(9)" <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur_cont=="Rather rural")),
                 
                 "(10)" <- lm(thermo.diff ~ resent_munis.std + id.std + cntry, data = subset(dat_full, dat_full$rur_cont=="Very rural")),
                 "(11)" <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + cntry, data = subset(dat_full, dat_full$rur_cont=="Very rural")),
                 "(12)" <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur_cont=="Very rural"))
)

modelsummary(TableA11, fmt = 3, stars = T,
             title = 'OLS regression results: Place-based affective polarisation on place-based resentment and place-based identity, by self-classified urban-rural residence (full variable).', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("resent_munis.std" = "Place-based resentment (Std.)",
                          "id.std" = "Place-based identity (Std.)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "migbckYes" = "Migration background (b.=no)",
                          "lr.std" = "Left-right (Std.)",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.19 Figure A12 ##########################################################
GALTAN <- lm(galtan_vote ~ thermo.diff*rur_cont + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = dat_full)

GALTAN <- predict_response(GALTAN, terms = c("thermo.diff", "rur_cont"))

FigA12_1 <- ggplot(GALTAN, aes(x = x, y = predicted, group = group)) +
  geom_line(aes(linetype = group)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = .2) +
  labs(y = "← GAL-TAN →",
       x = "Thermometer differential") +
  scale_fill_manual(values = c("black", "darkgrey", "grey60", "lightgrey", "grey30")) +
  scale_color_manual(values = c("black", "darkgrey", "grey60", "lightgrey", "grey30")) +
  scale_linetype_manual(values = c("solid", "dashed", "dotted", "dotdash","longdash")) +  
  theme_pubr() +
  theme(legend.position = "bottom", 
        panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title = element_blank(), 
        panel.spacing.x = unit(1.5, "lines"))

hist_theme <- theme_void() + theme(legend.position = "none")
dat_hist <- dat_full %>% dplyr::select(galtan_vote, thermo.diff, res, gender, age.std, educ, inc_dec, migbck, lr.std, cntry) %>% na.omit()

FigA12_2 <- ggplot(dat_hist, aes(x = thermo.diff, fill = res, linetype=res)) + 
  geom_density(alpha=.5) +
  scale_fill_manual(values = c("black", "darkgrey", "grey60", "lightgrey", "grey30")) +  
  hist_theme

png(file=here("FigureA12.png"),
    width = 8, height = 6, units = 'in', res = 800)
FigA12_2 + 
  FigA12_1 + plot_layout(byrow = FALSE, 
                         widths = c(1, 1), 
                         heights = unit(c(3, 1), c('cm', 'null')), ncol = 1)
dev.off()

# 6.20 Table A12 ############################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes"), ncol = 4)) 
attr(FE_add, "position") = 29

TableA12 <- list("(1)" <- lm(galtan_vote ~ thermo.diff.std*rur_cont + cntry, data = dat_full),
                 "(2)" <- lm(galtan_vote ~ thermo.diff.std*rur_cont + gender + age.std + educ + inc_dec + migbck + cntry, data = dat_full),
                 "(3)" <- lm(galtan_vote ~ thermo.diff.std*rur_cont + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = dat_full)
)

modelsummary(TableA12, fmt = 3, stars = T,
             title = 'OLS regression results: GAL-TAN voting on place-based affective polarisation, conditional on urban-rural self-classifications (full variable).', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("thermo.diff.std" = "Thermometer differential (Std.)",
                          "rur_contRather urban" = "Rather urban residence",
                          "rur_contRather rural" = "Rather rural residence",
                          "rur_contVery rural" = "Very rural residence",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "migbckYes" = "Migration background (b.=no)",
                          "lr.std" = "Left-right (Std.)",
                          "thermo.diff.std:rur_contRather urban" = "Thermometer differential (Std.) X Rather urban residence",
                          "thermo.diff.std:rur_contRather rural" = "Thermometer differential (Std.) X Rather rural residence",
                          "thermo.diff.std:rur_contVery rural" = "Thermometer differential (Std.) X Very rural residence",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.21 Figure A13 ##########################################################
inter1 <- interflex(estimator = 'binning', Y = "galtan_vote", D = "rur", X = "thermo.diff",
                    Z = c("gender", "age.std", "educ", "inc_dec", "migbck" ,"lr.std"),
                    FE = c("cntry"), 
                    data = dat_full, na.rm = TRUE,  base = "Urban residence") 

FigA13.1 <- plot(inter1, ylab = "Marginal effect of rural residence on \nGAL-TAN", xlab = "Thermometer differential", theme.bw = TRUE, Xdistr = "density")

inter2 <- interflex(estimator = 'binning', Y = "galtan_vote", D = "rur", X = "thermo.diff",
                    Z = c("gender", "age.std", "educ", "inc_dec", "migbck", "lr.std"),
                    FE = c("cntry"),
                    data = dat_full, na.rm = TRUE,  base = "Rural residence") 

FigA13.2 <- plot(inter2, ylab = "Marginal effect of urban residence on \nGAL-TAN", xlab = "Thermometer differential", theme.bw = TRUE, Xdistr = "density")

png(file=here("FigureA13.png"),
    width = 9, height = 6, units = 'in', res = 800)
ggarrange(FigA13.1, FigA13.2, ncol = 2, nrow = 1)
dev.off()

# 6.22 Figure A14 ##########################################################
thermo1 <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + immi.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence"))
summary(thermo1)
dat_thermo1 <- broom::tidy(thermo1)
dat_thermo1 %<>% filter(term=="resent_munis.std" | term=="id.std")
dat_thermo1$sample <- "Urban residence"

thermo2 <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + immi.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence"))
summary(thermo2)
dat_thermo2 <- broom::tidy(thermo2)
dat_thermo2 %<>% filter(term=="resent_munis.std" | term=="id.std")
dat_thermo2$sample <- "Rural residence"

dat_thermo_fin <- rbind(dat_thermo1, dat_thermo2)

dat_thermo_fin$term[dat_thermo_fin$term=="resent_munis.std"] <- "Place-based \nresentment (Std.)"
dat_thermo_fin$term[dat_thermo_fin$term=="id.std"] <- "Place-based \nidentity (Std.)"

dat_thermo_fin$sample <- as.factor(dat_thermo_fin$sample)
dat_thermo_fin$sample <- fct_rev(dat_thermo_fin$sample)
dat_thermo_fin$estimate <- as.numeric(dat_thermo_fin$estimate)
dat_thermo_fin$std.error <- as.numeric(dat_thermo_fin$std.error)
dat_thermo_fin$var <- dat_thermo_fin$term
dat_thermo_fin$var <- factor(as.factor(dat_thermo_fin$var), levels = c("Place-based \nidentity (Std.)", "Place-based \nresentment (Std.)"))


FigA14 <- ggplot(dat_thermo_fin, aes(x = var, y = estimate, ))
FigA14 <- FigA14 + geom_hline(yintercept = 0, color = "black", lty = 2)

FigA14 <- FigA14 + geom_linerange(aes(x = var, ymin = estimate - std.error*interval1,
                                      ymax = estimate + std.error*interval1, color = sample),
                                  lwd = 1, position = position_dodge(width = 1/2))  +
  scale_color_manual(name = "Subsample:", values = c("black", "red"))


FigA14 <- FigA14 + geom_pointrange(aes(x = var, y = estimate, ymin = estimate - std.error*interval2,
                                       ymax = estimate + std.error*interval2, color = sample, shape = sample),
                                   lwd = 1/2, position = position_dodge(width = 1/2)) +
  scale_color_manual(name = "Subsample:", values = c("black", "red")) +
  scale_shape_manual(name = "Subsample:", values = c(1, 2))


FigA14 <- FigA14 + xlab("Independent variable:") + ylab("Estimated coefficient \n Dependent variable: place-based affective polarisation") +  
  theme_pubr() +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_text(), panel.spacing.x = unit(1, "lines")) 

FigA14 <- FigA14 + coord_flip(ylim = c(-3.5, 10.5)) 

png(file=here("FigureA14.png"),
    width = 8, height = 3, units = 'in', res = 800)
FigA14 
dev.off()

# 6.23 Table A13 ###########################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes", "Yes"), ncol = 3)) 
attr(FE_add, "position") = 19

TableA13 <- list("(1)" <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + immi.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                 "(2)" <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + immi.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence"))
)

modelsummary(TableA13, fmt = 3, stars = T,
             title = 'OLS regression results: place-based affective polarisation on place-based resentment and place-based identity, by self-classified urban-rural residence (interactions).', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("resent_munis.std" = "Place-based resentment (Std.)",
                          "id.std" = "Place-based identity (Std.)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "migbckYes" = "Migrant background (b.=no)",
                          "immi.std" = "Immigration attitudes (Std.)",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.24 Figure A15 ##########################################################
GALTAN <- lm(galtan_vote ~ thermo.diff*rur + gender + age.std + educ + inc_dec + migbck + immi.std + cntry, data = dat_full)

GALTAN <- predict_response(GALTAN, terms = c("thermo.diff", "rur"))

FigA15_1 <- ggplot(GALTAN, aes(x = x, y = predicted, group = group)) +
  geom_line(aes(linetype=group)) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill=group), alpha = .2) +
  labs(y = "← GAL-TAN →",
       x = "Thermometer differential")  +
  scale_fill_manual(values = c("black", "grey")) + 
  theme_pubr() +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_blank(), panel.spacing.x = unit(1.5, "lines"))

hist_theme <- theme_void() + theme(legend.position = "none")
dat_hist <- dat_full %>% dplyr::select(galtan_vote, thermo.diff, rur, gender, age.std, educ, inc_dec, migbck, immi.std, cntry) %>% na.omit()

FigA15_2 <- ggplot(dat_hist, aes(x = thermo.diff, fill = rur, linetype=rur)) + 
  geom_density(alpha=.5) +
  scale_fill_manual(values = c("black", "grey")) +  
  hist_theme

png(file=here("FigureA15.png"),
    width = 8, height = 6, units = 'in', res = 800)
FigA15_2 + 
  FigA15_1 + plot_layout(byrow = FALSE, 
                         widths = c(1, 1), 
                         heights = unit(c(3, 1), c('cm', 'null')), ncol = 1)
dev.off()

# 6.25 Table A14 ###########################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes"), ncol = 2)) 
attr(FE_add, "position") = 21

TableA14 <- list("(1)" <- lm(galtan_vote ~ thermo.diff.std*rur + gender + age.std + educ + inc_dec + migbck + immi.std + cntry, data = dat_full)
)

modelsummary(TableA14, fmt = 3, stars = T,
             title = 'OLS regression results: GAL-TAN voting on place-based affective polarisation, conditional on urban-rural self-classifications (controlling for immigration attitudes).', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("thermo.diff.std" = "Thermometer differential (Std.)",
                          "rurRural residence" = "Rural residence (b.=urban residence)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "migbckYes" = "Migrant background (b.=no)",
                          "immi.std" = "Immigration attitudes (Std.)",
                          "thermo.diff.std:rurRural residence" = "Thermometer differential (Std.) X Rural residence",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.26 Figure A16 ##########################################################
AP_comp <- dat_full %>% dplyr::select(cntry, AP_Wag, thermo.in, thermo.out)

AP_comp$like_mean <- (AP_comp$thermo.in+AP_comp$thermo.out)/2
AP_comp$like_diff1 <- (AP_comp$thermo.in-AP_comp$like_mean)^2
AP_comp$like_diff2 <- (AP_comp$thermo.out-AP_comp$like_mean)^2

AP_comp$like_sum <- AP_comp$like_diff1+AP_comp$like_diff2
AP_comp$like_sum_div <- AP_comp$like_sum/2
AP_comp$AP_ur <- sqrt(AP_comp$like_sum_div)

AP_comp %<>% dplyr::select(cntry, AP_Wag, AP_ur)

grp.mean1 <- AP_comp %>%
  group_by(cntry) %>%
  dplyr::summarize(AP_Wag = mean(AP_Wag, na.rm=TRUE))

grp.mean2 <- AP_comp %>%
  group_by(cntry) %>%
  dplyr::summarize(AP_ur = mean(AP_ur, na.rm=TRUE))

grp.mean <- merge(grp.mean1, grp.mean2, by=c("cntry"))


AP_comp_lng <- melt(AP_comp, id.vars=c("cntry"))
AP_comp_lng$variable <- ifelse(AP_comp_lng$variable=="AP_Wag", "Affective partisan \npolarisation",
                               "Place-based \naffective polarisation")


grp.mean_lng <- melt(grp.mean, id.vars=c("cntry"))
grp.mean_lng$variable <- ifelse(grp.mean_lng$variable=="AP_Wag", "Affective partisan \npolarisation",
                                "Place-based \naffective polarisation")

png(file=here("FigureA16.png"),
    width = 9, height = 6, units = 'in', res = 800)
ggplot(AP_comp_lng, aes(x = value, fill = variable)) + 
  geom_density(alpha=.5) +
  scale_fill_manual(values = c("black", "grey")) + 
  labs(y = "Density",
       x = "Affective polarisation") + 
  facet_wrap(vars(cntry)) +
  geom_vline(data=grp.mean_lng, aes(xintercept=value, color = variable),
             linetype="dashed") + 
  scale_color_manual(values = c("grey", "black")) +
  theme_pubr() +
  theme(legend.position="bottom", panel.grid.major = element_line(color = "grey"),
        panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
        axis.line = element_line(colour = "black"),
        legend.title=element_blank(), panel.spacing.x = unit(1.5, "lines"))
dev.off()

AP_comp_corr <- AP_comp %>% select(AP_Wag, AP_ur)
cor(AP_comp_corr$AP_Wag, AP_comp_corr$AP_ur)

# 6.27 Table A15 ############################################################
t_tests <- list(
  CZ = t.test(thermo.diff ~ rur, data = dat_full %>% filter(cntry == "Czech Republic"), alternative = "less"),
  DK = t.test(thermo.diff ~ rur, data = dat_full %>% filter(cntry == "Denmark"), alternative = "less"),
  FR = t.test(thermo.diff ~ rur, data = dat_full %>% filter(cntry == "France"), alternative = "less"),
  DE = t.test(thermo.diff ~ rur, data = dat_full %>% filter(cntry == "Germany"), alternative = "less"),
  EL = t.test(thermo.diff ~ rur, data = dat_full %>% filter(cntry == "Greece"), alternative = "less"),
  HU = t.test(thermo.diff ~ rur, data = dat_full %>% filter(cntry == "Hungary"), alternative = "less"),
  IT = t.test(thermo.diff ~ rur, data = dat_full %>% filter(cntry == "Italy"), alternative = "less"),
  PL = t.test(thermo.diff ~ rur, data = dat_full %>% filter(cntry == "Poland"), alternative = "less"),
  ES = t.test(thermo.diff ~ rur, data = dat_full %>% filter(cntry == "Spain"), alternative = "less")
)

t_test_results <- data.frame(
  Country = names(t_tests),
  Estimate = sapply(t_tests, function(x) x$estimate[2] - x$estimate[1]),
  T_Statistic = sapply(t_tests, `[[`, "statistic"),
  Degrees_of_Freedom = sapply(t_tests, `[[`, "parameter"),
  P_Value = sapply(t_tests, function(x) format(x$p.value, digits = 4))
)

kable(t_test_results, "latex", booktabs = TRUE, digits = 3, caption = "Results of unpaired t-tests comparing thermometer differential scores 
between rural and urban residents per country (one-tailed).") %>%
  kable_styling(full_width = FALSE) 

# 6.27 Figure A17 ############################################################
bar_data <- as.data.frame.matrix(table(dat_full$res, dat_full$rur_cont))
bar_data <- bar_data %>%
  mutate(res = rownames(bar_data)) %>%
  pivot_longer(cols = -res, names_to = "rur", values_to = "count") %>%
  mutate(
    res = factor(res, levels = c(
      "A big city", 
      "The suburbs or outskirts of a big city", 
      "A town or small city", 
      "A country village", 
      "A farm or home in the countryside"
    )),
    rur = factor(rur, levels = c(
      "Very urban", 
      "Rather urban", 
      "Rather rural", 
      "Very rural"
    ), ordered = TRUE)  
  ) %>%
  group_by(rur) %>%
  mutate(percentage = count / sum(count) * 100)

png(file=here("FigureA17.png"),
    width = 9, height = 6, units = 'in', res = 800)
ggplot(bar_data, aes(x = rur, y = percentage, fill = res)) +
  geom_bar(stat = "identity", position = "stack", color = "black") +  
  scale_fill_brewer(
    palette = "Greys",
    name = "Current residence"  
  ) +
  labs(
    x = "Urban-rural (full variable)",
    y = "Percentage"
  ) +
  theme_pubr() +
  theme(
    legend.position = "bottom",
    legend.spacing.x = unit(0.5, "lines"),
    legend.text = element_text(size = 10),
    legend.key.width = unit(1, "cm"),
    legend.box = "horizontal",
    legend.box.margin = margin(t = 5),
    panel.grid.major = element_line(color = "grey"),
    panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
    axis.line = element_line(colour = "black"),
    panel.spacing.x = unit(1.5, "lines")
  ) +
  guides(fill = guide_legend(nrow = 2, byrow = TRUE)) +  
  scale_y_continuous(labels = scales::percent_format(scale = 1))
dev.off()

# 6.28 Table A16 ###########################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes", "Yes"), ncol = 3)) 
attr(FE_add, "position") = 7

TableA16 <- list("(1)" <- lm(thermo.diff ~ urban_move + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                 "(2)" <- lm(thermo.diff ~ rural_move + cntry, data = subset(dat_full, dat_full$rur=="Rural residence"))
)

modelsummary(TableA16, fmt = 3, stars = T,
             title = 'OLS regression results: place-based affective polarisation on movement indicators, by self-classified urban-rural residence.', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("urban_moveYes" = "Urban move (b.=no)",
                          "rural_moveYes" = "Rural move (b.=no)",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.29 Table A17 ###########################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes", "Yes"), ncol = 5)) 
attr(FE_add, "position") = 19

TableA17 <- list("(1)" <- lm(thermo.diff ~ urban_move + gender + age.std + educ + inc_dec + migbck + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                "(2)" <- lm(thermo.diff ~ urban_move + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                
                "(3)" <- lm(thermo.diff ~ rural_move + gender + age.std + educ + inc_dec + migbck + cntry, data = subset(dat_full, dat_full$rur=="Rural residence")),
                "(4)" <- lm(thermo.diff ~ rural_move + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence"))
)

modelsummary(TableA17, fmt = 3, stars = T,
             title = 'OLS regression results: place-based affective polarisation on movement indicators, by self-classified urban-rural residence (with controls).', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("urban_moveYes" = "Urban move (b.=no)",
                          "rural_moveYes" = "Rural move (b.=no)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "migbckYes" = "Migration background (b.=no)",
                          "lr.std" = "Left-right (Std.)",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.30 Table A18 ###########################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), ncol = 7)) 
attr(FE_add, "position") = 19

TableA18 <- list("(1)" <- lm(thermo.diff ~ resent_munis.std + id.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                "(2)" <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                "(3)" <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                
                "(4)" <- lm(thermo.diff ~ resent_munis.std + id.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence")),
                "(5)" <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + cntry, data = subset(dat_full, dat_full$rur=="Rural residence")),
                "(6)" <- lm(thermo.diff ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence"))
)

modelsummary(TableA18, fmt = 3, stars = T,
             title = 'OLS regression results: Place-based affective polarisation on place-based resentment and place-based identity, by self-classified urban-rural residence.', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("resent_munis.std" = "Place-based resentment (Std.)",
                          "id.std" = "Place-based identity (Std.)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "migbckYes" = "Migration background (b.=no)",
                          "lr.std" = "Left-right (Std.)",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.31 Table A19 ############################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), ncol = 7)) 
attr(FE_add, "position") = 25

TableA19 <- list("(1)" <- lm(thermo.diff ~ resent_munis.std*rur + id.std + cntry, data = dat_full),
                 "(2)" <- lm(thermo.diff ~ resent_munis.std*rur + id.std + gender + age.std + educ + inc_dec + migbck + cntry, data = dat_full),
                 "(3)" <- lm(thermo.diff ~ resent_munis.std*rur + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = dat_full),
                 
                 "(4)" <- lm(thermo.diff ~ resent_munis.std + id.std*rur + cntry, data = dat_full),
                 "(5)" <- lm(thermo.diff ~ resent_munis.std + id.std*rur + gender + age.std + educ + inc_dec + migbck + cntry, data = dat_full),
                 "(5)" <- lm(thermo.diff ~ resent_munis.std + id.std*rur + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = dat_full)
)


modelsummary(TableA19, fmt = 3, stars = T,
             title = 'OLS regression results: place-based affective polarisation on place-based resentment and place-based identity, by self-classified urban-rural residence (interactions).', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("resent_munis.std" = "Place-based resentment (Std.)",
                          "id.std" = "Place-based identity (Std.)",
                          "rurRural residence" = "Rural residence (b.=urban residence)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "migbckYes" = "Migration background (b.=no)",
                          "lr.std" = "Left-right (Std.)",
                          "resent_munis.std:rurRural residence" = "Place-based resentment (Std.) X Rural residence",
                          "id.std:rurRural residence" = "Place-based identity (Std.) X Rural residence",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.32 Table A20 ############################################################
TableA20_1 <- list("(1)" <- lm(thermo.diff ~ resent_munis.std + id.std, data = subset(dat_full, dat_full$cntry=="Czech Republic" & dat_full$rur=="Urban residence")),
                  "(2)" <- lm(thermo.diff ~ resent_munis.std + id.std, data = subset(dat_full, dat_full$cntry=="Czech Republic" & dat_full$rur=="Rural residence")),
                  "(3)" <- lm(thermo.diff ~ resent_munis.std + id.std, data = subset(dat_full, dat_full$cntry=="Denmark" & dat_full$rur=="Urban residence")),
                  "(4)" <- lm(thermo.diff ~ resent_munis.std + id.std, data = subset(dat_full, dat_full$cntry=="Denmark" & dat_full$rur=="Rural residence")),
                  "(5)" <- lm(thermo.diff ~ resent_munis.std + id.std, data = subset(dat_full, dat_full$cntry=="France" & dat_full$rur=="Urban residence")),
                  "(6)" <- lm(thermo.diff ~ resent_munis.std + id.std, data = subset(dat_full, dat_full$cntry=="France" & dat_full$rur=="Rural residence")),
                  "(7)" <- lm(thermo.diff ~ resent_munis.std + id.std, data = subset(dat_full, dat_full$cntry=="Germany" & dat_full$rur=="Urban residence")),
                  "(8)" <- lm(thermo.diff ~ resent_munis.std + id.std, data = subset(dat_full, dat_full$cntry=="Germany" & dat_full$rur=="Rural residence")),
                  "(9)" <- lm(thermo.diff ~ resent_munis.std + id.std, data = subset(dat_full, dat_full$cntry=="Greece" & dat_full$rur=="Urban residence")),
                  "(10)" <- lm(thermo.diff ~ resent_munis.std + id.std, data = subset(dat_full, dat_full$cntry=="Greece" & dat_full$rur=="Rural residence"))
                  
)

modelsummary(TableA20_1, fmt = 3, stars = T,
             title = 'OLS regression results: place-based affective polarisation on place-based resentment and place-based identity, by urban-rural residence (per country).', 
             coef_omit = "cntry", 
             coef_map = c("resent_munis.std" = "Place-based resentment (Std.)",
                          "id.std" = "Place-based identity (Std.)",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))


TableA20_2 <- list("(11)" <- lm(thermo.diff ~ resent_munis.std + id.std, data = subset(dat_full, dat_full$cntry=="Hungary" & dat_full$rur=="Urban residence")),
                  "(12)" <- lm(thermo.diff ~ resent_munis.std + id.std, data = subset(dat_full, dat_full$cntry=="Hungary" & dat_full$rur=="Rural residence")),
                  "(13)" <- lm(thermo.diff ~ resent_munis.std + id.std, data = subset(dat_full, dat_full$cntry=="Italy" & dat_full$rur=="Urban residence")),
                  "(14)" <- lm(thermo.diff ~ resent_munis.std + id.std, data = subset(dat_full, dat_full$cntry=="Italy" & dat_full$rur=="Rural residence")),
                  "(15)" <- lm(thermo.diff ~ resent_munis.std + id.std, data = subset(dat_full, dat_full$cntry=="Poland" & dat_full$rur=="Urban residence")),
                  "(16)" <- lm(thermo.diff ~ resent_munis.std + id.std, data = subset(dat_full, dat_full$cntry=="Poland" & dat_full$rur=="Rural residence")),
                  "(17)" <- lm(thermo.diff ~ resent_munis.std + id.std, data = subset(dat_full, dat_full$cntry=="Spain" & dat_full$rur=="Urban residence")),
                  "(18)" <- lm(thermo.diff ~ resent_munis.std + id.std, data = subset(dat_full, dat_full$cntry=="Spain" & dat_full$rur=="Rural residence"))
)

modelsummary(TableA20_2, fmt = 3, stars = T,
             title = 'OLS regression results: place-based affective polarisation on place-based resentment and place-based identity, by urban-rural residence (per country).', 
             coef_omit = "cntry", 
             coef_map = c("resent_munis.std" = "Place-based resentment (Std.)",
                          "id.std" = "Place-based identity (Std.)",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))


# 6.33 Table A21 ############################################################
TableA21_1 <- list("(1)" <- lm(thermo.diff ~ resent_munis.std*rur, data = subset(dat_full, dat_full$cntry=="Czech Republic")),
                  "(2)" <- lm(thermo.diff ~ id.std*rur, data = subset(dat_full, dat_full$cntry=="Czech Republic")),
                  "(3)" <- lm(thermo.diff ~ resent_munis.std*rur, data = subset(dat_full, dat_full$cntry=="Denmark")),
                  "(4)" <- lm(thermo.diff ~ id.std*rur, data = subset(dat_full, dat_full$cntry=="Denmark")),
                  "(5)" <- lm(thermo.diff ~ resent_munis.std*rur, data = subset(dat_full, dat_full$cntry=="France")),
                  "(6)" <- lm(thermo.diff ~ id.std*rur, data = subset(dat_full, dat_full$cntry=="France")),
                  "(7)" <- lm(thermo.diff ~ resent_munis.std*rur, data = subset(dat_full, dat_full$cntry=="Germany")),
                  "(8)" <- lm(thermo.diff ~ id.std*rur, data = subset(dat_full, dat_full$cntry=="Germany")),
                  "(9)" <- lm(thermo.diff ~ resent_munis.std*rur, data = subset(dat_full, dat_full$cntry=="Greece")),
                  "(10)" <- lm(thermo.diff ~ id.std*rur, data = subset(dat_full, dat_full$cntry=="Greece"))
                  
)

modelsummary(TableA21_1, fmt = 3, stars = T,
             title = 'OLS regression results: place-based affective polarisation on place-based resentment and place-based identity, by urban-rural residence (interactions; per country).', 
             coef_omit = "cntry", 
             coef_map = c("resent_munis.std" = "Place-based resentment (Std.)",
                          "id.std" = "Place-based identity (Std.)",
                          "rurRural residence" = "Rural residence (b.=urban residence)",
                          "resent_munis.std:rurRural residence" = "Place-based resentment (Std.) X Rural residence",
                          "id.std:rurRural residence" = "Place-based identity (Std.) X Rural residence",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))


TableA21_2 <- list("(11)" <- lm(thermo.diff ~ resent_munis.std*rur, data = subset(dat_full, dat_full$cntry=="Hungary")),
                  "(12)" <- lm(thermo.diff ~ id.std*rur, data = subset(dat_full, dat_full$cntry=="Hungary")),
                  "(13)" <- lm(thermo.diff ~ resent_munis.std*rur, data = subset(dat_full, dat_full$cntry=="Italy")),
                  "(14)" <- lm(thermo.diff ~ id.std*rur, data = subset(dat_full, dat_full$cntry=="Italy")),
                  "(15)" <- lm(thermo.diff ~ resent_munis.std*rur, data = subset(dat_full, dat_full$cntry=="Poland")),
                  "(16)" <- lm(thermo.diff ~ id.std*rur, data = subset(dat_full, dat_full$cntry=="Poland")),
                  "(17)" <- lm(thermo.diff ~ resent_munis.std*rur, data = subset(dat_full, dat_full$cntry=="Spain")),
                  "(18)" <- lm(thermo.diff ~ id.std*rur, data = subset(dat_full, dat_full$cntry=="Spain"))
)

modelsummary(TableA21_2, fmt = 3, stars = T,
             title = 'OLS regression results: place-based affective polarisation on place-based resentment and place-based identity, by urban-rural residence (interactions; per country).', 
             coef_omit = "cntry", 
             coef_map = c("resent_munis.std" = "Place-based resentment (Std.)",
                          "id.std" = "Place-based identity (Std.)",
                          "rurRural residence" = "Rural residence (b.=urban residence)",
                          "resent_munis.std:rurRural residence" = "Place-based resentment (Std.) X Rural residence",
                          "id.std:rurRural residence" = "Place-based identity (Std.) X Rural residence",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.34 Table A22 ############################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), ncol = 13)) 
attr(FE_add, "position") = 19

TableA22 <- list("(1)" <- lm(thermo.in ~ resent_munis.std + id.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                 "(2)" <- lm(thermo.in ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                 "(3)" <- lm(thermo.in ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                 
                 "(4)" <- lm(thermo.out ~ resent_munis.std + id.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                 "(5)" <- lm(thermo.out ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                 "(6)" <- lm(thermo.out ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
                 
                 "(7)" <- lm(thermo.in ~ resent_munis.std + id.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence")),
                 "(8)" <- lm(thermo.in ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + cntry, data = subset(dat_full, dat_full$rur=="Rural residence")),
                 "(9)" <- lm(thermo.in ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence")),
                 
                 "(10)" <- lm(thermo.out ~ resent_munis.std + id.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence")),
                 "(11)" <- lm(thermo.out ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + cntry, data = subset(dat_full, dat_full$rur=="Rural residence")),
                 "(12)" <- lm(thermo.out ~ resent_munis.std + id.std + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence"))
                 
)

modelsummary(TableA22, fmt = 3, stars = T,
             title = 'OLS regression results: in-group affect and out-group affect on place-based resentment and place-based identity, by self-classified urban-rural residence.', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("resent_munis.std" = "Place-based resentment (Std.)",
                          "id.std" = "Place-based identity (Std.)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "migbckYes" = "Migration background (b.=no)",
                          "lr.std" = "Left-right (Std.)",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.35 Table A23 ############################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes"), ncol = 4)) 
attr(FE_add, "position") = 21

TableA23 <- list("(1)" <- lm(galtan_vote ~ thermo.diff.std*rur + cntry, data = dat_full),
                "(2)" <- lm(galtan_vote ~ thermo.diff.std*rur + gender + age.std + educ + inc_dec + migbck + cntry, data = dat_full),
                "(3)" <- lm(galtan_vote ~ thermo.diff.std*rur + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = dat_full)
)

modelsummary(TableA23, fmt = 3, stars = T,
             title = 'OLS regression results: GAL-TAN voting on place-based affective polarisation, conditional on urban-rural self-classifications.', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("thermo.diff.std" = "Thermometer differential (Std.)",
                          "rurRural residence" = "Rural residence (b.=urban residence)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "migbckYes" = "Migration background (b.=no)",
                          "lr.std" = "Left-right (Std.)",
                          "thermo.diff.std:rurRural residence" = "Thermometer differential (Std.) X Rural residence",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))


GALTAN <- lm(galtan_vote ~ thermo.diff.std*rur + gender + age.std + educ + inc_dec +  migbck + lr.std + cntry, data = dat_full)
avg_slopes(GALTAN, variables = "thermo.diff.std", by = "rur")

# 6.36 Table A24 ############################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes"), ncol = 4)) 
attr(FE_add, "position") = 25

TableA24 <- list("(1)" <- lm(galtan_vote ~ thermo.diff.std*rur + resent_munis.std*rur + cntry, data = dat_full),
                 "(2)" <- lm(galtan_vote ~ thermo.diff.std*rur + resent_munis.std*rur + gender + age.std + educ + inc_dec + migbck + cntry, data = dat_full),
                 "(3)" <- lm(galtan_vote ~ thermo.diff.std*rur + resent_munis.std*rur + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = dat_full)
)

modelsummary(TableA24, fmt = 3, stars = T,
             title = 'OLS regression results: GAL-TAN voting on place-based affective polarisation, conditional on urban-rural self-classifications.', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("thermo.diff.std" = "Thermometer differential (Std.)",
                          "resent_munis.std" = "Place-based resentment (Std.)",
                          "rurRural residence" = "Rural residence (b.=urban residence)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "migbckYes" = "Migration background (b.=no)",
                          "lr.std" = "Left-right (Std.)",
                          "thermo.diff.std:rurRural residence" = "Thermometer differential (Std.) X Rural residence",
                          "rurRural residence:resent_munis.std" = "Place-based resentment (Std.) X Rural residence",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

GALTAN <- lm(galtan_vote ~ thermo.diff.std*rur + resent_munis.std*rur + gender + age.std + educ + inc_dec + lr.std + migbck + cntry, data = dat_full)
avg_slopes(GALTAN, variables = "thermo.diff.std", by = "rur")
avg_slopes(GALTAN, variables = "resent_munis.std", by = "rur")

# 6.37 Table A25 ############################################################
TableA25_1 <- list("(1)" <- lm(galtan_vote ~ thermo.diff.std*rur, data = subset(dat_full, dat_full$cntry=="Czech Republic")),
                   "(2)" <- lm(galtan_vote ~ thermo.diff.std*rur, data = subset(dat_full, dat_full$cntry=="Denmark")),
                   "(3)" <- lm(galtan_vote ~ thermo.diff.std*rur, data = subset(dat_full, dat_full$cntry=="France")),
                   "(4)" <- lm(galtan_vote ~ thermo.diff.std*rur, data = subset(dat_full, dat_full$cntry=="Germany")),
                   "(5)" <- lm(galtan_vote ~ thermo.diff.std*rur, data = subset(dat_full, dat_full$cntry=="Greece"))
)

modelsummary(TableA25_1, fmt = 3, stars = T,
             title = 'OLS regression results: GAL-TAN voting on place-based affective polarisation, conditional on urban-rural self-classifications (per country).', 
             coef_omit = "cntry", 
             coef_map = c("thermo.diff.std" = "Thermometer differential (Std.)",
                          "rurRural residence" = "Rural residence (b.=urban residence)",
                          "thermo.diff.std:rurRural residence" = "Thermometer differential (Std.) X Rural residence",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

TableA25_2 <- list("(6)" <- lm(galtan_vote ~ thermo.diff.std*rur, data = subset(dat_full, dat_full$cntry=="Hungary")),
                   "(7)" <- lm(galtan_vote ~ thermo.diff.std*rur, data = subset(dat_full, dat_full$cntry=="Italy")),
                   "(8)" <- lm(galtan_vote ~ thermo.diff.std*rur, data = subset(dat_full, dat_full$cntry=="Poland")),
                   "(9)" <- lm(galtan_vote ~ thermo.diff.std*rur, data = subset(dat_full, dat_full$cntry=="Spain"))
)

modelsummary(TableA25_2, fmt = 3, stars = T,
             title = 'OLS regression results: GAL-TAN voting on place-based affective polarisation, conditional on urban-rural self-classifications (per country).', 
             coef_omit = "cntry", 
             coef_map = c("thermo.diff.std" = "Thermometer differential (Std.)",
                          "rurRural residence" = "Rural residence (b.=urban residence)",
                          "thermo.diff.std:rurRural residence" = "Thermometer differential (Std.) X Rural residence",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.38 Table A26 ############################################################
TableA26_1 <- list("(1)" <- lm(galtan19 ~ thermo.diff.std*rur, data = subset(dat_full, dat_full$cntry=="Czech Republic")),
                   "(2)" <- lm(galtan19 ~ thermo.diff.std*rur, data = subset(dat_full, dat_full$cntry=="Denmark")),
                   "(3)" <- lm(galtan19 ~ thermo.diff.std*rur, data = subset(dat_full, dat_full$cntry=="France")),
                   "(4)" <- lm(galtan19 ~ thermo.diff.std*rur, data = subset(dat_full, dat_full$cntry=="Germany")),
                   "(5)" <- lm(galtan19 ~ thermo.diff.std*rur, data = subset(dat_full, dat_full$cntry=="Greece"))
)

modelsummary(TableA26_1, fmt = 3, stars = T,
             title = 'OLS regression results: GAL-TAN voting on place-based affective polarisation, conditional on urban-rural self-classifications (CHES 2019 data; per country).', 
             coef_omit = "cntry", 
             coef_map = c("thermo.diff.std" = "Thermometer differential (Std.)",
                          "rurRural residence" = "Rural residence (b.=urban residence)",
                          "thermo.diff.std:rurRural residence" = "Thermometer differential (Std.) X Rural residence",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

TableA26_2 <- list("(6)" <- lm(galtan19 ~ thermo.diff.std*rur, data = subset(dat_full, dat_full$cntry=="Hungary")),
                   "(7)" <- lm(galtan19 ~ thermo.diff.std*rur, data = subset(dat_full, dat_full$cntry=="Italy")),
                   "(8)" <- lm(galtan19 ~ thermo.diff.std*rur, data = subset(dat_full, dat_full$cntry=="Poland")),
                   "(9)" <- lm(galtan19 ~ thermo.diff.std*rur, data = subset(dat_full, dat_full$cntry=="Spain"))
)

modelsummary(TableA26_2, fmt = 3, stars = T,
             title = 'OLS regression results: GAL-TAN voting on place-based affective polarisation, conditional on urban-rural self-classifications (CHES 2019 data; per country).', 
             coef_omit = "cntry", 
             coef_map = c("thermo.diff.std" = "Thermometer differential (Std.)",
                          "rurRural residence" = "Rural residence (b.=urban residence)",
                          "thermo.diff.std:rurRural residence" = "Thermometer differential (Std.) X Rural residence",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.39 Table A27 ############################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes"), ncol = 4)) 
attr(FE_add, "position") = 23

TableA27 <- list("(1)" <- lm(galtan_vote ~ thermo.in.std*rur + thermo.out.std*rur + cntry, data = dat_full),
                 "(2)" <- lm(galtan_vote ~ thermo.in.std*rur + thermo.out.std*rur + gender + age.std + educ + inc_dec + migbck + cntry, data = dat_full),
                 "(3)" <- lm(galtan_vote ~ thermo.in.std*rur + thermo.out.std*rur + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = dat_full)
)


modelsummary(TableA27, fmt = 3, stars = T,
             title = 'OLS regression results: GAL-TAN voting on in-group affect and out-group affect, by self-classified urban-rural residence.', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("thermo.in.std" = "In-group affect (Std.)",
                          "thermo.out.std" = "Out-group affect (Std.)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "migbckYes" = "Migration background (b.=no)",
                          "lr.std" = "Left-right (Std.)",
                          "thermo.in.std:rurRural residence" = "In-group affect (Std.) X Rural residence",
                          "rurRural residence:thermo.out.std" = "Out-group affect (Std.) X Rural residence",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

GALTAN <- lm(galtan_vote ~ thermo.in.std*rur + thermo.out.std*rur + gender + age.std + educ + inc_dec + migbck + lr.std + cntry, data = dat_full)
avg_slopes(GALTAN, variables = "thermo.in.std", by = "rur")
avg_slopes(GALTAN, variables = "thermo.out.std", by = "rur")

# 6.40 Figure A18 ###########################################################
dat_ru_dens <- na.omit(dat_full[,c("popdens_n3","rur_cont")])

ur_dens <- dat_ru_dens %>% 
  group_by(rur_cont) %>% 
  summarise(
    mean_density = mean(popdens_n3, na.rm = TRUE),
    se = sd(popdens_n3, na.rm = TRUE) / sqrt(n()),  
    CI_lo = mean_density - interval1 * se,  
    CI_up = mean_density + interval1 * se   
  )

png(file=here("FigureA18.png"),
    width = 6, height = 3, units = 'in', res = 800)
ggplot(ur_dens, aes(x = rur_cont, y = mean_density)) +
  geom_linerange(aes(ymin = CI_lo, ymax = CI_up)) +
  geom_pointrange(aes(ymin = CI_lo, ymax = CI_up), size = 0.1) +
  labs(x = "Urban-rural (full variable)", y = "Average population density \nof NUTS-3 regions") +
  theme_pubr() +
  theme(
    legend.position = "bottom",
    legend.spacing.x = unit(0.5, "lines"),
    legend.text = element_text(size = 10),
    legend.key.width = unit(1, "cm"),
    legend.box = "horizontal",
    legend.box.margin = margin(t = 5),
    panel.grid.major = element_line(color = "grey"),
    panel.border = element_rect(color = "black", linewidth = 1, fill = NA),
    axis.line = element_line(colour = "black"),
    panel.spacing.x = unit(1.5, "lines")
  ) 
dev.off()

# 6.41 Table A28 ############################################################
avg_nuts2 <- dat_full %>%
  group_by(NUTS2) %>%
  summarise(count = n()) %>%
  summarise(mean_count = mean(count)) %>%
  pull(mean_count)
avg_nuts2

avg_nuts3 <- dat_full %>%
  group_by(NUTS3) %>%
  summarise(count = n()) %>%
  summarise(mean_count = mean(count)) %>%
  pull(mean_count)
avg_nuts3

dat_full$NUTS2 <- as.factor(dat_full$NUTS2)
dat_full$NUTS3 <- as.factor(dat_full$NUTS3)

empty1 <- lmer(thermo.diff ~ 1 + (1 | NUTS2), data = dat_full, REML = TRUE)
summary(empty1)
icc(empty1)

empty2 <- lmer(thermo.diff ~ 1 + (1 | NUTS3), data = dat_full, REML = TRUE)
summary(empty2)
icc(empty2)

empty3 <- lmer(galtan_vote ~ 1 + (1 | NUTS2), data = dat_full, REML = TRUE)
summary(empty3)
icc(empty3)

empty4 <- lmer(galtan_vote ~ 1 + (1 | NUTS3), data = dat_full, REML = TRUE)
summary(empty4)
icc(empty4)

# 6.42 Table A29 ############################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), ncol = 8)) 
attr(FE_add, "position") = 27

TableA29 <- list(
  "(1)" = lmer(galtan_vote ~ cntry +
                 (1 | NUTS2), data = dat_full, REML = TRUE),
  
  "(2)" = lmer(galtan_vote ~ thermo.diff.std*rur + cntry +
                 (1 | NUTS2), data = dat_full, REML = TRUE),
  
  "(3)" = lmer(galtan_vote ~ thermo.diff.std*rur + gender + age.std + educ + inc_dec + migbck + cntry + 
                 (1 | NUTS2), data = dat_full, REML = TRUE),
  
  "(4)" = lmer(galtan_vote ~ thermo.diff.std*rur + gender + age.std + educ + inc_dec + migbck + lr.std + cntry + 
                 (1 | NUTS2), data = dat_full, REML = TRUE),
  
  "(5)" = lmer(galtan_vote ~ thermo.diff.std*rur + gender + age.std + educ + inc_dec + migbck + lr.std + reg_unemp.std + cntry +
                 (1 | NUTS2), data = dat_full, REML = TRUE),
  
  "(6)" = lmer(galtan_vote ~ thermo.diff.std*rur + gender + age.std + educ + inc_dec + migbck + lr.std + reg_unemp.std + pop_delta.std + cntry + (1 | NUTS2), data = dat_full, REML = TRUE),
  
  "(7)" = lmer(galtan_vote ~ thermo.diff.std * rur + gender + age.std + educ + inc_dec + migbck + lr.std + reg_unemp.std + pop_delta.std + EQI + cntry + (1 | NUTS2), data = dat_full, REML = TRUE)
)

modelsummary(TableA29, fmt = 3, stars = T,
             title = '.', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c(
               "thermo.diff.std" = "Thermometer differential (Std.)",
               "rurRural residence" = "Rural residence (b.=urban residence)",
               "genderfemale" = "Gender (b.=male)",
               "age.std" = "Age (Std.)",
               "educHigh" = "Education (b.=low)",
               "inc_dec" = "Income (Deciles)",
               "migbckYes" = "Migration background (b.=no)",
               "lr.std" = "Left-right (Std.)",
               "thermo.diff.std:rurRural residence" = "Thermometer differential (Std.) X Rural residence",
               "reg_unemp.std" = "Regional unemployment (Std.)",
               "pop_delta.std" = "Population $\\Delta$ (Std.)",
               "EQI" = "EQI",
               "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))
